Lab 6
1. Expectation
Expectation does not exist: Cauchy Distribution
Make a sample from a Cauchy distribution and compute the sample mean and median.
z<- rcauchy(1000)
mean(z)## [1] 0.5208591
summary(z)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -198.2917 -0.8392 -0.0202 0.5209 0.9410 324.9145
y<- rnorm(1000,3,.5)
mean(y)## [1] 3.010735
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.250 2.652 3.019 3.011 3.348 4.420
Make many means and medians.
x.1 <- replicate(1000,mean(rcauchy(10000)))
x.2 <- replicate(1000, median(rcauchy(10000)))
summary(x.1) #summary of mean## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1925.7108 -0.8812 0.0439 -0.3628 1.0505 1953.1474
summary(x.2) #summary of median## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.0522713 -0.0107851 -0.0001419 0.0003810 0.0117740 0.0623587
Maybe the sample size is too small for computing the mean ?
Plot the ecdf of a suitable sample.
plot.ecdf(rcauchy(1000))Standard deviation of a Cauchy rv DOES NOT EXIST
boxplot(replicate(1000, sd(rcauchy(100))))2. Moments
Example1: Second moment of a binomial random variable
n = 30
p = .4
sum((0:n)^2*dbinom(0:n,n,p)) # E(X^2) = sum(x^2 * P(X=x))## [1] 151.2
mean(rbinom(10000,n,p)^2) # E(X^2)## [1] 151.0061
Example 2: Third central moment of a Poisson random variable
lambda = 4
n = 50
my.mean <- sum((0:n)*dpois(0:n, lambda)) # mu =E(X)=sum(x*P(X=x))
my.mean## [1] 4
sum(((0:n) - my.mean)^3*dpois(0:n, lambda)) #E((X-mu)^3) = sum((x-mu)^3 * P(X=x))## [1] 4
mean((rpois(100000,lambda) - my.mean)^3) #E((X-mu)^3)## [1] 3.80861
Example 3: Second central moment
- Standard deviation of a binomial distribution
n = 30
p = .4
sqrt(n*p*(1-p))## [1] 2.683282
sd(rbinom(10000,n,p))## [1] 2.687843
- Standard deviation of a normal distribution
my.mu = 1
my.sigma = 3.5
sd(rnorm(1000, my.mu, my.sigma))## [1] 3.465753
sqrt(mean((rnorm(1000,my.mu, my.sigma) - my.mu)^2)) #E((X-mu)^2))## [1] 3.522004
Example 4: Moments package
library(moments)
x <- rnorm(100)
## Compute the mean
moment(x) # E(X)## [1] 0.1765234
## Compute the 2nd centeral moment (= var)
moment(x, order=2,central=TRUE)## [1] 0.8559815
## Compute the 3rd absolute centeral moment
moment(x, order=3, central=TRUE, absolute=TRUE)## [1] 1.260265
3. Expected Value and Conditioning
- Expected value of \(X | X > 1/3\) if \(X \sim U(0,1)\)
x <- runif(10000)
mean(x[x>1/3]) #by simulation## [1] 0.6640886
2/3 # theorotical## [1] 0.6666667
- And the variance:
var(x[x>1/3])## [1] 0.03749393
1/27## [1] 0.03703704
- Conditional Distribution \(Y= X | X > 1/3\) if \(X \sim U(0,1)\)
x <- runif(1000)
y <- x[x>(1/3)]plot.ecdf(x)plot.ecdf(y)Example 5: Mixed Distributions
Let \(X.0 \sim U(0,3)\), \(X.1 \sim Poisson(4)\), and let \(X = X.0 + X.1\)
Make an empirical cdf of \(X\)
Make an empirical cdf of \(X|X < 8\) in the same plot
Make an empirical cdf of \(X.1|X<8\)
N = 10000
X.0 <- runif(N,0,3)
X.1 <- rpois(N,4)
X <- X.0 + X.1plot.ecdf(X)
X.2 <- X[X<8]
mycdf = ecdf(X.2)
t <- seq(0,10,by = .01)
lines(t,mycdf(t), col = 2)plot.ecdf(X.1[X<8])4. Simulation of bivariate random variables
Example 6: Joint Distribution
\(X1 > 0\), \(X2 > 0\), \(X1 + X2 < 1\), uniform density = 2
Also plot it
N = 10000 # number of initial simulation
set.seed(101)
mydf.0 <- data.frame(X1 = runif(N), X2 = runif(N))
mydf.0 <- mydf.0[rowSums(mydf.0) < 1,]
head(mydf.0)## X1 X2
## 2 0.04382482 0.55931317
## 5 0.24985572 0.37856509
## 8 0.33346714 0.52437539
## 12 0.70687474 0.02394005
## 15 0.45512059 0.06166735
## 19 0.41166683 0.15480003
plot(mydf.0, pch = 46)Example 7: to illustrate Law of Total Probability
- First simulate large number of conditions \(X \in A_j\). These occur as in a geometric distribution with p = 0.5.
_Because geometric p.m.f \(=p(1-p)^{x-1}=(1/2)(1/2)^{j-1}=(1/2)^j\)
n = 100000
A.1 <- rgeom(n,.5)
table(A.1)## A.1
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 49857 25097 12570 6243 3165 1507 800 378 202 78 48 28 12
## 13 14 15 16 17
## 6 2 2 3 2
- For each value j of this condition, simulate a uniform random variable in the interval \((j,j+1)\). Then make a histogram.
X.5 <- runif(n, min = A.1, max = A.1+1)
hist(X.5, prob = T, breaks =0:20)- Compute the sample mean:
mean(X.5)## [1] 1.499157
3/2## [1] 1.5
Example 8: Numerical Example
Suppose \((X,Y)\) have the joint pmf
\[ f_{XY}(1,1)= .1, \quad f_{XY}(1,2) = .2, \quad f_{XY}(2,1) = .1, \quad f_{XY}(2,2) = .6, \]
Find the expectations of \(X\) and \(Y\) \[ E(X) = 1.7, \quad E(Y) = 1.8 \]
Find the conditional distributions of \(X\), conditioned on the possible values \(y\) of \(Y\):
\[ \begin{aligned} f_{X|Y}(1|Y=1) &= .5, \quad f_{X|Y}(2|Y=1) = .5\\ f_{X|Y}(1|Y=2) &= .25, \quad f_{X|Y}(2|Y=2) = .75\\ \end{aligned} \]
__Conditional expectations_
\[ E(X|Y=1) = 1.5, E(X|Y=2) = 1.75 \] Formula for \(E(X|Y)\)
\[ E(X|Y) = \begin{cases} 1.5 \quad \text{if} \; Y=1 \\ 1.75 \quad \text{if} \; Y=2 \end{cases} \]
5. Marginal distribution, continuous case.
Example 9:
- Simulation of \((U,V)\)
mydf.0 = data.frame(U = runif(10000), V = runif(10000))
plot(mydf.0,cex = .1)- Simulation of \((X,Y)\)
x = sapply(1:10000, function(j){min(mydf.0[j,])})
y = sapply(1:10000, function(j){max(mydf.0[j,])})
# Can also be done with apply()
mydf.1 = data.frame(x =x, y = y)
plot(mydf.1,pch = 1,cex = .1)- For the marginal distribution of X, just make a probability histogram
hist(mydf.1$x, prob = T)
abline(a = 2, b = -2, col = 2) It shows the straight-line shape that is predicted by theory, $ 2-2x$
- marginal distribution of Y
hist(mydf.1$y, prob = T)
abline(a = 0, b = 2, col = 2) Past Exam Question: Joint distributions and Conditional Expectations.
A pdf is defined by
\[ \begin{aligned} f(x,y) &= \begin{cases} C(x+2y) \quad (0<y<1, 0<x<2)\\ 0 \quad (otherwise) \end{cases} \end{aligned} \]
Find the value of \(C\).
Find the Marginal density of \(X\)
Find the Marginal density of \(Y\)
Find the conditional density function of of \(Y|X=x\)
Find the Conditional expectation \(E(Y|X)\)